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ASSUMED-STRESS HYBRID ELEMENTS WITH DRILLING DEGREES OF FREEDOM 
FOR NONLINEAR ANALYSIS OF COMPOSITE STRUCTURES 
- FINAL REPORT FOR NAG-1 -1505 - 


Norman F. Knight, Jr. 

Department of Aerospace Engineering 
Old Dominion University 
Norfolk, VA 23529-0247 

SUMMARY 

The goal of this research project is to develop assumed-stress hybrid elements with rota- 
tional degrees of freedom for analyzing composite structures. During the first year of the 
three-year activity, the effort was directed to further assess the AQ4 shell element and its 
extensions to buckling and free vibration problems. In addition, the development of a com- 
patible 2-node beam element was to be accomplished. The extensions and new develop- 
ments were implemented in the Computational Structural Mechanics Testbed COMET. An 
assessment was performed to verify the implementation and to assess the performance 
of these elements in terms of accuracy. During the second and third years, extensions to 
geometrically nonlinear problems were developed and tested. This effort involved working 
with the nonlinear solution strategy as well as the nonlinear formulation for the elements. 

This research has resulted in the development and implementation of two additional ele- 
ment processors (ES22 for the beam element and ES24 for the shell elements) in COMET. 
The software was developed using a SUN workstation and has been ported to the NASA 
Langley Convex named blackbird. Both element processors are now part of the baseline 
version of COMET. 

The initial grant (NAG-1 -1374) had a period of performance of January 10, 1992 to Janu- 
ary 14, 1993 with an approved no-cost extension until May 15, 1993 while the principal in- 
vestigator was affiliated with Clemson University. Dr. Alexander Tessler was the Contract- 
ing Officer’s Technical Representative (COTR). In January 1 993, the principal investigator 
moved to Old Dominion University and a new grant activity was established 
(NAG-1-1505). This grant had a period of performance of April 10, 1993 until May 15, 
1996 with an approved no-cost extension until December 31, 1996. Initially, Dr. Tessler 
remained the COTR for NAG-1 -1 505; then Mr. Keith Norwood assumed that role. When 
Mr. Norwood left NASA Langley, Dr. W. Jefferson Stroud became the COTR fora short peri- 
od and then Dr. Tessler re-assumed that role. 



OVERVIEW OF GRANT 


NAG-1-1374 (Clemson First Year: 1/10/92 to 5/15/93) 


Grant Objectives 

Grant Accomplishments 

Assess AQ4 formulation, implementation, capabilities 

Developed symbolic approach for deriving these ele- 
ments based on MAPLE; compared computational ef- 
forts (1 journal paper); 1 MS on beam element; 

Develop compatible beam element 

Develop compatible triangular element 

Deferred; told to emphasize beam element first 

Demonstrate combined use of elements 

Sobel used elements on Grumman panel 

Extend to buckling problems 

Developed geometric stiffness matrices 

Extend to free vibration problems 

Developed consistent mass matrices 

Extend to geometrically nonlinear problems 

Deferred 

Utilize GEP in COMET 

ES22 (beam element) and ES24 (quad element) 




NAG-1-1505 (ODU First Year: 4/16/93 to 5/16/94) 


Grant Objectives 

Grant Accomplishments 

Explore alternative stress fields 

1 MS on quad element; 1 journal paper on formulation 

Derive diagonal mass coefficients 

HZR mass lumping for quad and beam 

Assess alternate stress recovery procedures 

Stress parameters; differentiating displacement field 

Redevelop beam element for GCP formulation 

Force/moment resultants in LAUB & GCP different 

Develop further test cases for quad element 

Skewed plate buckling; cylinder buckling & vibration; 
pear-shaped cylinder 

Develop compatible triangular element 

Deferred 


Demonstrate combined use of elements 


Supplement for ADR Work (Summer 1993) 

Explore use of ADR on parallel computers 

Effective algorithm developed; 4 journal papers, 1 PhD; 
presentation to CSB in Aug. 1993 

Supplement for Interface Element Work (Spring Semester 1994) 

Extend to elements with drilling freedoms 

Software changes 

Extend to curved 1-D interfaces 

1 SDM paper 

Assess solvers 

Software access 



NAG-1-1505 (ODU Second Year: 5/16/94 to 12/31/95) 


Grant Objectives 

Grant Accomplishments 

Review of NL_STATIC_1 in preparation for nonlinear 
analysis work 

Tutorial given to CSB on 1/6/95; problems found with 
ES6 

Extend beam to geometrically nonlinear problems 

Develop internal force vector; Demonstrate planar and 
spatial response prediction; 1 MS; 1 SDM paper; placed 
2nd at Mid-Atlantic AIAA student paper competition in 
April 1995; ES22 documentation 

Extend quad to geometrically nonlinear problems 

Develop internal force vector; Demonstrated accuracy 
of low-order corotational approach for problems; 

Develop compatible triangular element 

Initiated but then deferred due to difficulties encoun- 
tered with quad element for Raasch problem and with 
nonlinear problems having high number of halfwaves. 

Explore use of interface element to surfaces 

Work transferred; separate grant with Amin pour as PI 



















































BACKGROUND 


Finite element technology continues to be a major research topic in shell analysis studies. 
Although the finite element method is well established in the engineering community, sev- 
eral areas of element research are active and new contributions are being made. These 
research areas are focused on the element formulation, eliminating element sensitivity to 
mesh distortion, and providing a compatible set of element for modeling aerospace struc- 
tures. The first element technology research area is associated with the element formula- 
tion itself. Generally the formulation is denoted as either C° or C 1 . Formulations denoted 
as C° include the effects of transverse shear deformation which implies that the out-of- 
plane deflections and the bending rotations are approximated independently. In a C° ele- 
ment, the out-of-plane deflections and the bending rotations may be approximated using 
the same interpolation or shape functions (e. g., bilinear or biquadratic functions). Formula- 
tions denoted as C 1 are based on classical formulations which assume zero transverse 
shear flexibility. In a C 1 element, the out-of-plane deflections and the bending rotations are 
generally coupled in that a cubic interpolation is assumed for the out-of-plane deflections 
(e. g., Hermitian shape functions). In general, C° elements are more readily implemented 
than C 1 elements; however, C° elements are also more susceptible to rank deficiencies 
which may cause element “locking”. While reduced integration (either uniform or selective) 
has been used to alleviate element “locking,” it may also introduce spurious modes in the 
finite element analysis solution. Element formulations also differ in their starting point from 
the weak-form of the variational statement. For example, formulations based on the princi- 
ple of minimum potential energy results in displacement-based elements. However, 
through the use of alternate variational statements (e.g., Hellinger-Reissner principle), as- 
sumed stress hybrid formulations can be developed. In a hybrid formulation, the stress field 
across the element domain is assumed in such a way that the equilibrium equations are 
satisfied on the interior of the element and the displacement field is assumed only along 
the element boundaries. Elements based on a hybrid formulation are generally more accu- 
rate than displacement-based elements; however, they do pose difficulties in extending 
these elements to nonlinear problems. 

The second element technology research area is associated with eliminating mesh distor- 
tion sensitivity. Distorted meshes, where the element shape deviates significantly from a 
planar square, are common occurrences in most structural analysis problems. These 
meshes are caused by the complexity of the structural geometry, by the use of mesh transi- 
tion schemes to reduce the size of the computational problem, and by including nonlinear 
effects which may transform a regular mesh into a distorted mesh due to large deflections 
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and rotations. The effect of mesh distortion on shell response predictions has been re- 
ported by Stanley et al. (1986), Park and Stanley (1986), and Stanley (1989). Recently new 
formulations have been emerging which reduce or eliminate the effects of mesh distortion 
from the element formulation. One such formulation is the assumed natural-coordinate 
strain or ANS formulation developed by Park and Stanley (1986). Their 9-node ANS quad- 
rilateral shell element has been shown to be nearly insensitive to mesh distortion. However, 
this element fails the patch test. Another formulation is the 4-node assumed-stress hybrid 
shell element (AQ4) developed by Aminpour (1990a,b) which uses the same bending dis- 
placement field approximations as Tessler and Hughes (1 983). The AQ4 element has rota- 
tional degrees of freedom that are normal to the element surface or drilling freedoms which 
are used to improve the inplane displacement field approximations. The AQ4 element has 
been demonstrated to be accurate, pass both the membrane and bending patch tests, is 
nearly insensitive to mesh distortion, does not lock, and is invariant with respect to node 
numbering. 

The third element technology research area is associated with providing a compatible set 
of elements for modeling aerospace structures. Developers of robust shell elements must 
also be concerned with how the new elements will be used by the structural analysis com- 
munity. Large-scale complexity applications require a compatible set, or family, of ele- 
ments for modeling and analysis. This family of elements should be capable of represent- 
ing both 1 -D and 2-D structures and should provide quadrilateral as well as triangular shell 
elements along with compatible beam elements. Having these elements available, ana- 
lysts are provided the modeling flexibility needed to represent aerospace structures. 

ORIGINAL RESEARCH OBJECTIVES 

The specific research objectives initially proposed for this grant were as followed: 

• During the first year (9/1/91 to 8/31/92), an independent assessment of the per- 
formance of the AQ4 element will be made using test cases involving composite 
and built-up structures. In addition a 3-node triangular shell element (AQ3) and 
a 2-node beam element (AQ2) will be developed and implemented in COMET 
using the GEP. These elements will be compatible with the AQ4 shell element. 

• During the second year (9/1/92 to 8/31/93), the formulation of each of these ele- 
ments will be extended to address buckling problems. Specific test cases will 
be posed to verify the implementation and the also the robustness of the ele- 
ments. In addition, the modeling flexibility provided by this family of elements will 
be demonstrated (i.e., analyze a model which incorporates the AQ2, AQ3, and 
AQ4 elements. 

• During the third year (9/1/93 to 8/31/94), the formulation of each of these ele- 
ments will be extended to address geometrically nonlinear analyses. Specific 
test cases will be posed to verify the implementation and the also the robustness 
of the elements. 
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However, these objectives were modified early in the first year after discussions with Dr. 
Jerrold Housner during the first grant review on May 28, 1992. The specific objectives for 
the first year were to concentrate the effort on the 4-node quadrilateral shell element and 
the 2-node beam element and provide capabilities for linear stress, buckling and free vibra- 
tion problems for both elements. 

FIRST-YEAR PROGRESS REPORT 

The progress made during the first year is divided into two parts. One part deals with the 
further assessment of the AQ4 shell element and its extensions to linear buckling and free 
vibration problems. This research was conducted in part by Mr. Govind Rengarajan and 
resulted in his Master’s thesis (Rengarajan, 1 993). The second part deals with the develop- 
ment of an AQ4-compatible 2-node beam element with capabilities for linear stress, buck- 
ling and free vibration analyses. This research was conducted in part by Mr. Venkateshwar 
Deshpande and resulted in his Master’s thesis (Deshpande, 1993). 

Shell Element Research 

The work related to the 4-node shell element began with Aminpour’s AQ4 element 
(1990a,b) and implemented in COMET (Stewart, 1989) as element processor ES8/AQ4 
using the generic element processor (Stanley and Nour-Omid, 1990). A family of as- 
sumed-stress hybrid shell elements was installed as element processor ES24 in COMET. 
Three areas of research were investigated: (i) use of symbolic computations to accelerate 
the element computations; (ii) extensions to linear buckling and free vibration analyses; 
and (iii) use of bubble functions and higher order stress approximations. The basic results 
for each area will now be summarized. 

Use of Symbolic Computations - Assumed-stress hybrid elements are usually 
associated with higher computational cost due to the fact that a matrix inverse is needed 
in the evaluation of certain element matrices. For example the element stiffness matrix for 
the AQ4 shell element is a 24x24 matrix that is computed as followed: 


[^e] 24x24 = [^j r [/Y e ] “ 1 [7^] 
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where 


[He] 22x22 


L^eJ 22x24 



[P] T [D][P]dA 


[P] T [L][N]dA 


The matrix [D] is the compliance matrix from the strain-stress relations; the matrix [L] is 
the linear differential operator matrix from the strain-displacement relations; the matrix [P] 
is the matrix of stress field approximations; and the matrix [A(| is the matrix of displacement 
field approximations or shape functions. The need to invert the matrix [H e ] for each element 
and to evaluate two integrals over the area of the element significantly increases the ele- 
ment computation time. Traditionally these operations are performed numerically. That 
is, Gaussian quadrature is used to evaluate each entry in the matrices [ H e ] and [T e ], and 
then an explicit matrix inverse of [H e ] is computed. In this research, symbolic computations 
using the MAPLE software system (Char et al., 1991) are used to evaluate the two integrals 
exactly. Explicit symbolic expressions for [H e ] inverse could not be obtained even when 
represented as the adjoint and determinant of [H e ]. However, since the explicit inverse is 
not required but rather the matrix product [H e ]- ! [3^], this product is readily determined in 
an efficient manner by solving a system of equations with multiple right-hand-side vectors. 
Results obtained to-date indicate that the symbolic implementation of AQ4 (called 
ES24/A4S1) executes approximately twice as fast as the numerically integrated original 
version (ES8/AQ4). 


Extensions to Linear Buckling and Free Vibration -The original formulation of the AQ4 
shell element did not address the problems of linear buckling or free vibration. To address 
these problems, the formulation is extended and required the development of the geomet- 
ric stiffness matrix and the consistent mass matrix. The geometric stiffness matrix is based 
on the membrane stress resultant prebuckling stress state and the nonlinear Green-La- 
grange membrane strains. The geometric stiffness matrix is evaluated using Gauss quad- 
rature after the integrand of each entry is evaluated symbolically using MAPLE. Numerical 
integration is required for the geometric stiffness matrix for the case of general quadrilateral 
elements wherein the Jacobian is a function of the natural coordinates of the element. As 
such, the integrand involves ratios of polynomials and cannot be easily integrated symboli- 
cally. The mass matrix is a consistent mass matrix and includes the effects of rotary inertia. 
The integrand of the mass matrix involves only polynomials of the natural coordinates and 
as such can be easily integrated symbolically using MAPLE. 
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Performance of this element for problems involving linear buckling and free vibration are 
reported by Rengarajan (1993). Comparison with other 4-node shell elements in COMET 
(4_ANS, 4_HYB, 4_STG) indicate that the AQ4 shell element formulation is substantially 
better both for buckling and vibration of flat and curved structures. 

Beam Element Research 

The work related to the beam element is to develop an AQ4-compatible 2-node beam ele- 
ment based on the Hellinger-Reissner variational principle. A 2-node assumed-stress hy- 
brid 3-D beam element was installed as element processor ES22 in COMET using the ge- 
neric element processor. Three areas of research were investigated: (i) development of 
the beam element; (ii) use of symbolic computations to accelerate the element computa- 
tions; and (iii) extensions to linear buckling and free vibration analyses. The basic results 
for each area will now be summarized. 

Development of Beam Element - Based on the Hellinger-Reissner variational principle, 
a 2-node assumed-stress hybrid 3-D beam was developed. The element has six degrees 
of freedom per node (3 translations and 3 rotations) plus six independent stress parame- 
ters. The effects of transverse shear deformation and rotary inertia are included in the for- 
mulation; however, the effects of cross-sectional warping are neglected. The beam ele- 
ment is implemented as element processor ES22 in COMET. Details of the formulation, 
implementation and testing of this element is reported by Deshpande (1 993). Performance 
of this beam element for linear stress analysis is compared to that of the beam element in 
processor ES6/E21 0 - namely the STAGS 21 0 C 1 2-node beam element. In general, the 
two elements give comparable results for thin slender beams whether straight or curved. 
However; as the span-to-thickness ratio decreases, the new element provides a signifi- 
cant improvement in solution since the effects of transverse shear deformation are included 
in the formulation. 

Use of Symbolic Computations - Similar to the shell elements, assumed-stress hybrid 
beam elements are usually associated with higher computational cost due to the fact that 
a matrix inverse is needed in the evaluation of certain element matrices. For example, the 
element stiffness matrix for the ES22/ABS2 beam element is a 12x12 matrix that is com- 
puted as followed: 

fc] 12X12 = iTe] T [He]-'[T e ] 
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where 


[tf,] 6 x6 


[^6x12 



[P] T [D][P]dx 


[P] T [L][N]dx 


The matrix [£>] is the compliance matrix from the strain-stress relations; the matrix [L] is 
the linear differential operator matrix from the strain-displacement relations; the matrix [/>] 
is the matrix of stress field approximations; and the matrix [A/] is the matrix of displacement 
field approximations or shape functions. The need to invert the matrix [H e ] for each element 
and to evaluate two integrals along the length of the element significantly increases the ele- 
ment computation time. In this research, symbolic computations using the MAPLE soft- 
ware system are used to evaluate the two integrals exactly. Explicit symbolic expressions 
for [H e ] inverse were obtained in terms of the adjoint and determinant of [H e ] since this ma- 
trix is only a 6x6 matrix. This complete symbolic form of the element is implement as ele- 
ment ABS2 in element processor ES22. In order to assess, a version of the beam element 
was developed (ES22/ABN2) which evaluates the integrals symbolically and determines 
the matrix product [He]- 1 [7;] by solving a system of equations with multiple right-hand-side 
vectors. It was determined that the numerical evaluation of the matrix product is approxi- 
mately 5% faster than the approach using a symbolic form of the matrix inverse and then 
a matrix-matrix multiply. This is due in part to the complexity of the FORTRAN expressions 
for the symbolic forms of the adjoint and determinant of [ H e ]. 


For the 4-node shell element, the symbolic evaluation of [H e ]~ } is not possible because 
the symbolic expressions become unwieldy for a 24x24 matrix. However, since the explicit 
inverse is not required but rather the matrix product [H e ]- ! [7^], this product is readily deter- 
mined in an efficient manner by solving a system of equations with multiple right-hand-side 
vectors. The “cost” of using this advanced multi-field element has been assessed and 
compared with other 4— node elements implemented in COMET. Significant computational 
savings have been accrued by using a combination of symbolic integration procedures and 
numerical analysis routines. These results are reported in a paper published in Commu- 
nications in Numerical Methods for Engineering. 


Extensions to Linear Buckling and Free Vibration - The formulation of the beam ele- 
ment was extended to linear buckling and free vibration analyses and required the develop- 
ment of the geometric stiffness matrix and the consistent mass matrix. The geometric stiff- 
ness matrix is based on the axial force resultant prebuckling stress state and the nonlinear 
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Green-Lagrange axial strain. The geometric stiffness matrix is evaluated symbolically us- 
ing MAPLE. The mass matrix is a consistent mass matrix and includes the effects of rotary 
inertia. The integrand of the mass matrix is also easily integrated symbolically using MA- 
PLE. 

Performance of this element for problems involving linear buckling and free vibration are 
reported by Deshpande (1 993). Comparison with the ES6/E210 2-node beam element in 
COMET indicate that the ES22/ABN2 beam element formulation does require more ele- 
ments for convergence of higher modes and that convergence is to a lower value since 
transverse shear deformation and rotary inertia effects are included in the ES22/ABN2 for- 
mulation. Also the span-to-thickness ratio decreases, the new element provides a signifi- 
cant improvement in solution accuracy for linear buckling and free vibration problems. 

SECOND-YEAR PROGRESS REPORT 

This grant represents a continuation of the research effort initiated by the principal investi- 
gator while at Clemson University underNASA Grant No. NAG-1-1374. Afterthe principal 
investigator relocated to the Old Dominion University, the research effort was renewed and 
initiated at ODU under NASA Grant No. NAG— 1—1505 in April 1993. As such, this research 
thrust has been underway for approximately two years (one year at Clemson University and 
one year at ODU). Progress under this phase of the grant was reported in an oral presenta- 
tion to the Computational Mechanics Branch at NASA Langley Research Center on June 
30, 1994. Copies of the presentation material were given to NASA at that time. 

The research effort has focussed primarily on completing a more thorough assessment of 
the assumed stress hybrid formulation by examining different assumed stress fields, on 
deriving diagonal mass coefficients for the beam and quadrilateral shell elements, and for 
preparing to move into the area of geometric nonlinear response prediction. Specific re- 
search objectives for the second year (first year at ODU) were as followed: 

• Extend the 2-node assumed-stress hybrid beam element to handle geometri- 
cally nonlinear problems using both low-order and high-order corotation ap- 
proaches. 

• Validate the combined use of the assumed-stress hybrid 1 -D and 2-D elements 
for modeling built-up structures. 

• Extend the 2-D assumed-stress hybrid shell elements to handle geometrically 
nonlinear problems using both low— order and high— order corotation ap- 
proaches. 

• Develop a compatible 3-node assumed-stress hybrid shell element with drilling 
degrees of freedom with capabilities for linear stress, buckling and vibration 
analyses. 
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• Examine the drilling freedom formulations based on Allman— type shape func- 
tions and on the unsymmetric stress tensor. 

Element Research 

The first activity was to re-formulate the beam and shell elements to use the GCP or Gener- 
ic Constitutive Processor for the constitutive matrix. Previously, the processor LAU was 
used for the shell element, and the processor LAUB was used for the beam element. The 
GCP uses the same format as processor LAU for shell elements, but a different convention 
than LAUB for the beam element. This required a complete re-derivation of the beam ele- 
ment kernel routines followed by an extensive period of testing and validation. At present, 
the shell element can handle any composite laminate and accommodates shear correction 
factors as implemented by the GCP. The beam element expects the GCP to provide the 
constitutive stiffness terms typically called EA, El, and GJ. The beam element processor 
can handle any type of beam (isotropic or laminated) provided that the stiffness terms can 
be defined using the GCP. 

The second activity was to develop and implement a diagonal mass matrix for the beam 
element. This was performed using the Hinton-Rock-Zienkiewicz approach. Options are 
now provided for either a consistent mass matrix or a lumped diagonal mass matrix. Free 
vibration test cases were performed and reported in Mr. W. Scott Carron’s master’s thesis 
(Carron, 1995). 

The third activity was to investigate the stress recovery procedure for these elements. 
Since these elements are assumed stress-hybrid elements, the stress parameters can be 
computed and saved for each element or recovered from the element nodal displace- 
ments. As such, the element stresses (actually stress resultants) can be computed any- 
where in the element accurately without differentiating the element nodal displacement 
field. However, it does require either storing or re-computing certain element terms during 
the solution process. As an alternative approach, a stress recovery procedure was investi- 
gated in which the element strains are computed by differentiating the nodal shape func- 
tions for the displacements - basically discarding the assumed-stress functions. The nu- 
merical results reported in Carron’s thesis (Carron, 1995) indicate that the stress recovery 
procedure based on the assumed stress field was always better than that obtained by dif- 
ferentiating the nodal displacements. This came as no surprise, but served as a good ex- 
ample of a side benefit of these elements (improved displacement solutions and accurate 
stress recovery). 
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Table 1 . Family of assumed-stress hybrid shell elements. 


Element 

Name 

Additional 

Modifications 

Number of 
displacement 
dof 

Number of 
Stress Parameters 

Inplane 

Out-of- 

Plane 

Membrane 

Bending 

A4S1 

Symbolic Version 
of AQ4 

12 

12 

9 

13 

A4S2 

Modified Mem- 
brane Stress Field 

12 

12 

11 

13 

A4S3 

Bubble functions 
for out-of-plane 
only 

12 

15 

9 

17 

A4S4 

Bubble functions 
for inplane only 

14 

12 

11 

13 

A4S5 

Bubble functions 
for both 

14 

15 

11 

17 

A4S6 

Bubble functions 
for out-of-plane 
only 

12 

15 

9 

19 


Use of Bubble Functions 

The finite element approximations for the 4-node element can be expected to improve by 
adding degrees of freedom at the center of the element. The contribution of these degrees 
of freedom should not affect the edges of the elements. One such approximation is called 
a “bubble function” which has a value of unity at the element center in the natural coordi- 
nate frame and zero along all edges of the element. Addition of the bubble function for the 
translational and bending rotation degrees of freedom are considered independently and 
require modification of the stress approximations. These modifications led to the develop- 
ment of a family of assumed-stress hybrid shell elements with drilling degrees of freedom 
denoted as A4S/' where / ranges from 1 to 6. The features of these elements are summa- 
rized in Table 1 . 

The performance of each element of this family was assessed for linear stress analysis 
only. The results indicate that overall the A4S1 element is very robust and accurate. Ele- 
ments A4S2 and A4S4 are also quite good but increase the element-level computations 
and did not show any improvement over the performance of the A4S1 element. The re- 
maining elements (A4S3, A4S5, and A4S6) appear to be quite stiff and sensitive to extreme 
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Table 2. Performance of elements using the MacNeal— Harder cantilever beam. 


A : 

Extension, B : In-plane Shear, C 

: Out-of-plane Shear, D : Twist 

Load 

4_ANS 

4J4SC 

4.STG 

A4S1 1 A4S2 

A4S3 | A4S4 | A4S5 

A4S6 

j~~ Rectangular-Shag 

jed Elements 

A 

BBI 

0.995 


tjfcKKf 

0.998 

0.998 

0.988 

0.988 

0.988 

B 


0.904* 

■ 


0.993 

0.993 

0.993 

0.993 

0.993 

C 


0.986 

■ I 


0.981 

0.981 

0.981 

0.981 

0.981 

D 

0.856 

0.941 

0.680 

1.009 

1.009 

1.009 

1.009 

1.009 

0.858 

Trapezoid 

al-Shaped Elements 

A 


ffiTjTjTjH 



0.998 


0.998 

0.998 

0.998 

B 


rJjtS 

■ 

gjj^ 3 

0.985 


0.986 

0.986 

0.986 

C 

WtM 


t 

ISmSI 

0.969 

0.968 

0.969 

0.968 

0.961 

D 

0.843 

0.951 

t 

1.007 

1.007 

1.004 

1.007 

1.004 

0.856 


Parallelogram-Shaped Elements 


A 

0.966 

0.996 

9.989 


0.998 

0.998 

0.998 

0.998 

0.998 

B 

0.324 

0.080* 

0.794 

0.977 

0.972 

0.977 

0.977 

0.977 

0.977 

C 

0.939 

0.977 

0.991 

0.980 

0.980 

0.980 

0.980 

0.980 

0.979 

D 

0.798 

0.945 

0.677 

1.007 

1.007 

0.999 

1.007 

0.999 

0.846 


* Q4S results for these cases axe 0.993, 0.988, and 0.986, 
respectively. 

f Produces a singular stiffness matrix. 

mesh distortion. The assessment of different assumed stress fields and the performance 
of the element for linear stress, buckling, and free vibration analyses has been completed 
and is reported in a paper in the International Journal for Numerical Methods in Engineer- 
ing. As a result, future work is directed only towards the A4S1 element formulation. 


Drilling Freedoms 

In first-order shear deformation plate theory, the kinematic relations involve three transla- 
tions and two bending rotations - the rotation normal to the plate surface, the so-called 
drilling rotation does not enter the kinematic relations. Hence at an element level, the dril- 
ling degree of freedom is not present; however, at the global or assembled level, the de- 
grees of freedom associated with the structure include the drilling rotation. While assem- 
bling the element stiffness matrices to form the global stiffness matrix, if the elements 
connected to a particular node happen to be coplanar, then the drilling rotation of the node 
is not resisted and the global stiffness matrix is singular. This singularity could be avoided 
if the drilling degree of freedom is included at the element level. Generally the drilling de- 
grees of freedom are suppressed at the beginning of an analysis since they do not enter 
the kinematic description of the problem. There have been numerous attempts in formulat- 
ing elements with drilling degrees of freedom. Two main approaches have been success- 
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ful. One approach is to introduce the drilling rotations in the functional as independent vari- 
ables, while the second approach is to introduce the drilling rotations in the finite element 
displacement approximations. 

In the first approach, the true drilling or inplane rotation is defined in terms of the skew-sym- 
metric part of the displacement gradient as 



Reissner (1965) presented a mixed variational principle in which this definition of the in- 
plane rotation was relaxed and then later enforced in the variational statement as a 
constraint through the use of Lagrange multipliers. On imposing the stationary conditions 
on the functional, the Lagrange multipliers were identified as the skew— symmetric part of 
the stress tensor. Reissner’s formulation introduced the drilling rotation independent of the 
inplane displacements. Hughes and Brezzi (1989) have shown that Reissner’s mixed for- 
mulation holds good in a continuous case, but becomes unstable if discrete approximations 
of the continuous space is made using standard interpolation functions. They modified Re- 
issner’s functional to include an additional term which stabilizes the function in discrete 
approximations. Ibrahimbegovic, Taylor, and Wilson (1990) used the Hughes and Brezzi 
formulation to develop a quadrilateral membrane element using independent interpolation 
fields for rotations and Allman-type shape functions for the inplane displacements. Later 
Ibrahimbegovic and Wilson (1991) developed triangular and quadrilateral flat shell ele- 
ments along the same lines. The drilling rotations which occur in Allman-type shape func- 
tions are not true drilling rotations rather they are more like rotQtionsI conn&ctors (Amin- 
pour, 1990). Hence lura and Atluri (1992) argue that in using Allman— type shape functions 
compatibility will not be satisfied directly along nodal lines of the shell element. They have 
developed an element where the rotations are interpolated from the true drilling rotations 
evaluated from the skew-symmetric part of the displacement gradient at the nodes. These 
issues are still unresolved. 

In the second approach, the drilling rotation was introduced through the inplane displace- 
ment field approximations (e.g., Allman, 1984; Bergan and Felippa, 1985). Allman (1984) 
used a quadratic function for the normal component of the displacement and a linear func- 
tion for the tangential component of the displacement along each edge of the element. 
Cook (1986) noted that the same shape functions can be derived by eliminating the transla- 
tional displacements at the midpoint of an edge in favor of the comer drilling rotations. He 
applied it to an 8-node quadrilateral membrane element to obtain a 4-node quadrilateral 
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element. Soon other investigators followed suit including the work of Aminpour (1990) 
which was based on an assumed— stress hybrid formulation. These elements have been 
shown to be quite robust for linear stress analysis problems — particularly for membrane 
problems. The present work involves extensions of the work of Aminpour (1 990) to assess 
linear buckling, free vibration and geometrical nonlinearities. 

Application Studies 

Skewed Laminated Plate - Buckling of a skewed laminated plate subjected to combined 
axial compression and inplane shear was studied to consider the performance of the A4S1 
(new AQ4) element for buckling analysis of a composite structure with skewed elements. 
The plate is skewed by 45°, the laminate is a [±45/90/0] s laminate, and a 31x31 mesh of 
nodes is used to compute the reference solution. This mesh and a contour plot of the trans- 
verse deflection component of the buckling mode shape are shown in Figure 1 . Conver- 
gence of the buckling coefficient as a function of the number of nodes along each edge of 
the plate is shown in Figure 2 for two sets of boundary conditions (i.e., all edges clamped 
and all edges simply support). In both cases, the present A4S1 element converges quickly 
and correlates well with the solution obtained using the STAGS shell element named 
4_STG. 

Pear-Shaped Cylinder - Linear stress and buckling analyses were performed on the 
pear-shaped cylinder shown in Figure 3. It has a uniform thickness and is simply sup- 
ported. Only one— quarter of the cylinder is modeled, and the loading is uniform end short- 
ening. This noncircular cylindrical shell problem poses certain difficulties due to the 
changes in curvature around the shell. Comparison of the transverse deflection normalized 
by the shell thickness as a function of the circumferential coordinate is shown in Figure 4. 
Results are obtained using both the 4ANS element and the A4S1 (new AQ4) element. 
These results are in very good agreement. Linear buckling analyses were also carried out 
on this cylinder, and good correlation with reference solution was obtained. The buckling 
modes for the lowest two modes are shown in Figure 5. The buckling loads are normalized 
by the buckling loads obtained using the 9ANS element and a 7x51 mesh (i.e., 4.393 and 
6.483, respectively). 
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Figure 1. Refined finite element mesh and buckling mode shape contour plot. 
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Figure 3. Pear-shaped cylinder: geometry and material properties. 



Figure 4. Comparison of transverse deflection distribution around the shell. 


Mode 1 


(P cr) converged 

A4S1 =1.013 
4ANS = 1 .021 



Mode 2 


(Per) converged 

A4S1 =1.016 
4ANS = 1 .032 



Figure 5. Comparison of linear buckling analysis results. 
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Circular Cylindrical Shell - To assess the eigenvalue analysis capability of the AQ4 shell 
element, linear buckling and vibration analyses were carried out on an elastic isotropic thin 
circular cylinder. The cylinder is simply supported at both ends and has an R/t ratio of 1 00. 
The geometric stiffness matrix and the consistent mass matrix used in these computations 
are formed symbolically and then implemented into COMET using the GEP . In the analysis, 
a one— eighth symmetric model is used; however, for visualizing the eigenmode shapes, 
a full cylinder is considered. 

Comparison of the buckling results using the 4-node shell elements in COMET is given 
Rengarajan (1 993). The buckling mode shape is shown in Figure 6(a). These finite element 
results are compared with the theoretical solution and indicate that the ES24/A4S1 or AQ4 
element produces significantly better results than the other elements for the same discreti- 
zation due in part to the assumed-stress field features of the element. For the coarse 
mesh, the use of the STAGS 41 0 element (ES5/E41 0) resulted in triggering a “mechanism” 
in the mode shape. 

Comparison of the vibration results using the 4-node shell elements is given Rengarajan 
(1993). These results indicate that the AQ4 element produces better results than the other 
shear— flexible elements due in part to the higher— order assumed— displacement field by in- 
cluding the drilling freedoms. The AQ4 assumed-displacement field is still of lower order 
than the STAGS element which neglects transverse shear effects. The vibration mode 
shape corresponding to the lowest vibration frequency is shown in Figure 6(b). 




(a) Buckling mode shape 


(b) Vibration Mode Shape. 


Figure 6. Eigenmode shapes. 
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Grumman Shear Panel — Another assessment problem was associated with a large finite 
element mode' of a shear buckling problem of interest to Grumman Aerospace Corpora- 
tion. A significant amount of time was allocated to applying the beam and quadrilateral shell 
elements to this NASA-customer application problem - Grumman shear buckling compos- 
ite panel. The original analysis was carried out by Dr. Larry Sobel of Grumman using the 
Grumman version of STAGS and then using the STAGS elements in COMET. Difficulties 
were encountered by Dr. Sobel in performing the COMET analysis. These new COMET 
results did not correlate with his existing STAGS results, and it was unclear whether the 
difficulties were associated with the STAGS-element implementation in COMET, the non- 
linear analysis procedure, or some other modeling issues. Eventually, an implementation 
error associated with the sign of the twisting moment was detected in the STAGS shell ele- 
ment kernel routines implemented in COMET as processor ES5. 

However, when the Grumman composite shear panel was analyzed using the elements 
developed under this research grant, linear stress and buckling results were obtained that 
matched both other analytical results as well as results obtained from experiment. This 
complex model involved a combination of the 2— node beam and 4— node shell elements in 
a single finite element model and demonstrate their combined use in large finite element 
simulations. 

Nonlinear Analysis Procedure 

Next steps were taken to validate the nonlinear solution procedure available in COMET’S 
CLAMP procedure nl_static_l. Solution to small test problems appeared to be readily 
obtained and accurate. However, application to NASA finite element models of large stiff- 
ened composite panels led to convergence problems in the postbuckled region of the re- 
sponse. Numerous test cases and example problems were analyzed; yet the same conver- 
gence difficulties could not be duplicated except with the large model. This large model had 
many complicating features including multi— point constraints. A definitive resolution to this 
problem was not found; however, issues related to the modeling strategy were raised and 
were to be considered by the NASA engineers. Even though this activity consumed a large 
amount of time and produced little direct benefit for the research thrust of the grant, it did 
provide a thorough assessment of the nonlinear solution procedure and its performance 
which is an integral part of the research. The nonlinear analysis procedure has now been 
thoroughly tested and several improvements made. A tutorial presentation on the nonlin- 
ear solution procedure was given to personnel from the Computational Mechanics Branch 
on January 6, 1995 following a grant review presentation that same day. 
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The basic strategy for extending the assumed-stress hybrid beam and shell elements to 
handle geometric nonlinearities is through the use of the Generic Element Processor (or 
GEP) of COMET. The GEP provides a variety of utilities including both those for constitutive 
modeling (namely, the Generic Constitutive Processor or GCP) and those for large arbitrary 
rotations (namely, the corotation utilities). By using the corotational utilities, these elements 
will be extended to problems involving geometric nonlinearities by two approaches. The 
first approach is the low-order corotation approach that involves linear strain-displace- 
ment relations at the element level and the large displacement/rotation capability at the 
global level through the corotation utilities. This is accomplished by formulating and imple- 
menting the internal force vector for this family of elements in addition to the capabilities 
provided under NASA Grant No. NAG— 1—1374. The low— order corotation approach is 
straightforward for the beam element since the drilling degrees of freedom are actual de- 
grees of freedom that enter the formulation through the kinematics of the problem. Howev- 
er, this approach has not been demonstrated for plate and shell elements with drilling de- 
grees of freedom. 

The second approach is the high— order corotation approach that involves nonlinear strain- 
displacement relations at the element level and the large displacement/rotation capability 
at the global level through the corotation utilities. In this case, the internal force vector and 
material stiffness matrix will need to be formulated in order to account for the nonlinear 
strain— displacement relations and implemented. The high— order corotational approach 
and the total Lagrangian approach for assumed— stress hybrid elements poses some chal- 
lenges that will increase the computational effort at the element level — similar to that of dis- 
placement-based elements. 

Adaptive Dynamic Relaxation Supplement 

A supplement to the grant enabled an exploratory study of the adaptive dynamic relaxation 
algorithm for three— dimension hyperelastic nonlinear structures on massively parallel pro- 
cessing (MPP) systems. This method exploits the advantages of explicit time integration 
for both static and transient dynamic analyses of structures. The ADR algorithm is de- 
scribed in a three-part paper published in the international journal for Computer Methods 
in Applied Mechanics and Engineering. The algorithm has been shown to be very scalable 
as the number of processors increases, has minimal memory requirements for each pro- 
cessor, and minimizes the interprocessor communication. For example, on the Intel Touch- 
stone Delta MPP system with 512 processors, a relative speed-up of 386 was obtained. 
Results from this study are summarized in a paper published in the International Journal 
for Numerical Methods in Engineering. 
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Interface Element Supplement 

Upon the hiring of Dr. Mohammad A. Aminpour by ODU in January 1 994, a supplement task 
was initiated under this grant related to interface element technology. This task allowed 
Dr. Aminpour to begin his research upon arrival at ODU. Dr. Aminpour reported on this as- 
pect of the work independently. 

THIRD-YEAR PROGRESS REPORT 

This grant represents a continuation of the research effort initiated by the principal investi- 
gator at Old Dominion University under NASA Grant No. NAG-1-1 505 which was renewed 
in May 1994. The activity continued until May 1995 and then was extended to December 
31 , 1995 under a no-cost extension. As such, this research thrust has been underway for 
approximately three years (one year at Clemson University and two years at ODU). Prog- 
ress under this phase of the grant was reported in two oral presentations to the Computa- 
tional Mechanics Branch at NASA Langley Research Center on January 6, 1995 and on 
December 12,1 995. Copies of the presentation material were given to NASA at the time 
of presentation. 

Specific research objectives for the third year (second year at ODU) were as followed: 

• Extend the 2-node assumed-stress hybrid beam element to handle geometri- 
cally nonlinear problems using both low-order and high-order corotation ap- 
proaches. 

• Extend the 2-D assumed-stress hybrid shell elements to handle geometrically 
nonlinear problems using both low-order and high-order corotation ap- 
proaches. 

• Validate the combined use of the assumed-stress hybrid 1-D and 2-D elements 
for modeling built-up structures. 

• Develop a compatible 3-node assumed-stress hybrid shell element with drilling 
degrees of freedom with capabilities for linear stress, buckling and vibration 
analyses. 


Element Research 

The element research performed during this period focussed on the extensions needed for 
geometrically nonlinear analysis for both the shell and beam elements. The approach tak- 
en was to exploit the corotational utilities and and existing element subroutines. To do this, 
only the internal force vector needed to be computed and used along with the linear materi- 
al stiffness matrix and the geometric stiffness matrix. The internal force vector is computed 
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as 



This form is readily available since the matrix [T] is computed for the linear material stiffness 
matrix and the vector of unknown stress parameters {0} are computed at the element level. 

In the preceding approach, the local element strain-displacement relations involve only the 
linear terms and the the nonlinear behavior is treated through the corotational formulation. 
In the corotational formulation, a local element coordinate system is established and ele- 
ment deformations are measure relative to the local system which corotates with the ele- 
ment through the global coordinate system. This approach is called low-order corotation. 
This approach works well provided that local element rotations are small. Typically results 
can be improved by refining the mesh which essentially adds more corotating coordinate 
systems. 

As an alternative, the local element strain-displacement relations could include the nonlin- 
ear terms as well as the linear terms. In so doing, each element would be able to treat large 
local rotations than in the low-order corotational formulation before refining the mesh 
would be necessary. However, this high-order corotational formulation for the assumed- 
stress hybrid elements requires the solution of a nonlinear system at the element level 
which significant increases the computational cost of the element. Presently, an assess- 
ment of the low-order corotational formulation for nonlinear analysis using the new beam 
and shell elements indicates that implementing the nonlinear strain-displacement rela- 
tions at the element level is not needed. 


Application Studies 

Many applications have been analyzed during this phase of the research. The applications 
studies involving the beam element are given by Carron (1995). Only two beam problems 
are considered herein: 45-degree circular bend problem and the clamped-hinged deep 
circular arch. The application studies involving the shell element are given here for se- 
lected applications although many more have been studied to insure the reliability of the 
element and its implementation. Three shell problems are considered: the elastica prob- 
lem, the hinged-free cylindrical panel, and the Raasch challenge problem. The final ap- 
plication problem involved a model of the HSCT using the AQ4 element for its 4-node 
quadrilateral elements. 
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Beam Cross-Section 



R=100.0 in. 
E=10 7 psi 
v=0.0 


8 2-node elements 


Figure 7. 45-Degree Circular Bend: Geometry and Properties. 



Figure 8. Deformed Shapes and Results Comparisons. 


45-degree Circular Bend Problem - The 45-degree circular bend problem shown in Fig- 
ure 7 is a severe test of the corotational approach and the 2-node beam element because 
of the three-dimensional aspects of the deformation process. The beam is originally in the 
X-Y plane and is pulled in the Z direction. The finite element model had only eight 2-node 
beam elements along the beam’s length, and the low-order corotational approach is used. 
The results are shown graphically on Figure 8 for two values of load and compared with 
the solutions obtained by Bathe and Bolourchi (1979). The position of Point A in the XYZ 
coordinate system from Bathe and Bolourchi (1 979) is compared with the position obtained 
using the present 2-node element, and the results are quite good. Comparison to other 
investigators is given by Carron (1995). 
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R=100 in. 
t=l. in. 
b=l . in. 
E=12xl0 6 psi 
v=0.0 


Figure 9. Clamped-hinged circular arch. 

Clamped-hinged Circular Arch problem - The circular arches considered herein have 
the same geometry as those studied by DaDeppo and Schmidt (1975) and shown in Figure 
9. The following notation is used: P is the downward concentrated point load acting at the 
crown of the arch at Point C; 2a is the subtending angle of the arch; R is the centroidal radius 
of curvature; E is the modulus of elasticity; v is Poisson’s ratio; I is the moment of inertia; 
U is the horizontal deflection of the crown point C; and V is the vertical deflection of the 
crown point C. The rectangular cross-sectional dimensions and elastic modulus are se- 
lected such that El = 1 0 6 lb-in. 2 Two sets of boundary conditions are considered. One set 
is asymmetric and are different at the two ends of the arch (one is clamped; the other is 
hinged), while the other set is symmetric with both ends being clamped. In both cases, the 
applied load Pat the arch crown (at Point C) is normalized such that p = pr 2 /ei . The most 
common geometry considered by other researchers is the case for a subtending angle 2a 
of 21 5° and asymmetric boundary conditions. This geometry corresponds to a deep circu- 
lar arch. The solutions obtained by DaDeppo and Schmidt (1975) are based on Euler’s 
nonlinear theory of an inextensible elastica with no restrictions on the magnitude of the 
deflections. For this subtending angle, the normalized collapse load F„from DaDeppo and 
Schmidt (1975) is 8.97. Numerical solutions are obtained using the present element with 
the low-order corotational formulation and the modified Newton-Raphson solution proce- 
dure with arc-length control. 

The load versus deflection curves obtained for the 20-element model are shown in Figure 
10 and compared with the analytical solutions of DaDeppo and Schmidt (1975). The 
normalized applied load P is shown as a function of the vertical deflection \/and the horizon- 
tal deflection U at Point C normalized by the centroidal radius of the arch R. These results 
clearly indicate the accuracy of the present element and the solution approach. Also shown 
in Figure 8 are the results obtained from a linear analysis which indicate the unconservative 
nature of such a prediction for nonlinear softening systems. 
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Figure 10. Load vs. deflections for clamped-hinged deep circular arch (2a=215°). 

Using the 20-element model, the subtending angle 2a is varied to determine its effect on 
the structural behavior of such elastic circular arches. Results of this parametric study 
along with the analytical solutions obtained by DaDeppo and Schmidt (1 975) are presented 
in Figure 1 1 . The solutions obtained using the present element are in good agreement with 
the reference solutions. As the subtending angle increases, a cusp in the collapse load 
versus subtending angle curve shown in Figure 11 occurs at approximately 2a=210°. On 
either side of this cusp, a smooth curve appears to describe the structural behavior. 

Examining the deformed geometries for arches with different subtending angles 2a on ei- 
ther side of the cusp reveals different structural behavior. For 2a=120°, the arch exhibits 
a prominent local bending behavior near the arch crown. The arch does not wrap around 
the hinged support and yet continues to snap through. For 2a=260°, the arch exhibits mi- 
nor local bending nearthe crown. However, the arch does exhibit a strong tendency to wrap 
around the hinged support in an asymmetric sidesway mode. For2a=215°, the arch exhib- 
its behavior that is more complex. As the load increases from zero, the arch deforms under 
local bending for a load level well below the collapse load. Just prior to reaching the col- 
lapse load, the asymmetry of the arch boundary conditions, as well as a local bending re- 
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60 80 100 120 140 160 180 200 220 240 260 

Subtending Angle 2a, degrees 

Figure 11 . Collapse load versus subtending angle for clamped-hinged circular arch. 


sponse or buckle near the crown, begin to be amplified. Near the collapse load, the arch 
begins to wrap around the hinged support. 

Several observations were made based on these studies reported by Carron (1 995). First, 
as the value of the subtending angle increases, the slope of these load versus deflection 
curves prior to collapse tends to increase. Also as the angle increases, the value of the 
collapse load decreases except for values near the cusp in Figure 11 . Near the cusp, the 
load-deflection curve exhibits a rapid increase in stiffness as the collapse load is ap- 
proached. For arches with a subtending angle greater than 210°, the nonlinear results indi- 
cate that, after very large pre-collapse deflections (when the crown has deflected vertically 
to the level of the supports), no further increase in the horizontal deflection U at the crown 
occurs. 
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10 in. 



M = 2n in.-lb/in. 


E=120 psi; v=0.3; thickness, t=1.0 in. 



Figure 12. Elastica Problem: geometry, material properties, and loading. 


Elastica Problem - The elastica problem is solved as another nonlinear problem with large 
deflections and rotations. The geometry, material properties and loading are shown in Fig- 
ure 12. The nonlinear elastic response of this cantilevered beam bent into a circle by an 
applied end moment as shown in Figure 12 is characterized by the moment versus end 
deflection response shown in Figure 13. The response is described in terms of the applied 
moment and the end deflection, a vertical tangent is approached as the end deflection 
approaches 1 .2 times the beam length. Three solutions are presented in Figure 11 . The 
reference solution is taken as the solution obtained using a 5x1 mesh of 9-node ANS ele- 
ments (ES1/EX97). Solutions were also obtained using 4-node shell elements with a 1 0x2 
mesh using the 4-node ANS element and the AQ4 element. All results are essentially the 
same through the loading range considered. 
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X Ref. solution (11x3 9ANS) 
▼11x3 AQ4 (ES24/A4S1 ) 

• 11x3 4ANS, NLGEOM=2 


Moment. 

M/2ti 



Figure 13. Elastica Problem with an End Moment: Moment vs. End Deflection. 


Free b 



R=2540 mm 
L=254 mm 

thickness, t=6.35 mm 
R/t=400 

E=3. 10275 kN/mm 2 
v=0.3 

<|>=0.1 radians 


Figure 14. Hinged-free cylindrical panel: geometry, material properties and loading. 
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Hinged— Free Cylindrical Panel — The hinged— free cylindrical panel subjected to a center 
transverse point load is considered to validate the nonlinear aspects of the AQ4 shell ele- 
ment based on the low-order corotational formulation. This structure is shown in Figure 
14. The structural response exhibits a limit point, a snap— through behavior, and a snap- 
back behavior. 

The nonlinear response of the panel is shown in Figure 15 for a midpoint on the free edge 
(Point b). Two meshes are considered, and the reference solution is taken to be that ob- 
tained using a 5x5 mesh of 9ANS elements. The results in Figure 15(a) indicate the re- 
sponse obtained using two meshes of the AQ4 element developed as part of this research. 
For both meshes, nearly the same nonlinear behavior is obtained. The only difference ap- 
pears to be after the shell’s free edge has snapped through to an inverted position. The 
results in Figure 1 5(b) indicate the response obtained using a 5x5 mesh of 4ANS elements 
with both the low— order and high-order corotational formulation (i.e., NLGEOM=1 and 
NLGEOM=2, respectively). These results appear to bracket the reference solution. 

The nonlinear response of the panel is also shown in Figure 16 for the center point (Point 
a). Two meshes are considered, and the reference solution is taken to be that obtained us- 
ing a 5x5 mesh of 9ANS elements.The results in Figure 16(a) indicate the response ob- 
tained using two meshes of the AQ4 element. For both meshes, nearly the same nonlinear 
behavior is obtained. Differences appear after the limit point as the snap— back response 
occurs. In this region of the nonlinear response, severe local changes in curvature occur 
and the 5x5 mesh of 4-node elements with the low-order corotation approach is not suffi- 
cient to reproduce the reference solution. Adding more elements (i.e., essentially adding 
more local reference frames) has the effect of refining the solution to obtain the reference 
solution. As an alternative, using the high-order corotational approach should have the 
same effect. The results in Figure 1 6(b) indicate the response obtained using a 5x5 mesh 
of 4ANS elements with both the low-order and high-order corotational formulation. These 
results are similar to the reference solution up until snap-back occurs. The two solutions 
exhibit a different character than anticipate from previous results. 
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X Ref. solution (5x5 9ANS) 
w 5x5 AQ4 (ES24/A4S1) 

• 7x7 AQ4 (ES24/A4S1) 



Load, 

kN 
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0.50 
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0.00 
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Transverse Deflection at Point b, mm 


(a) AQ4 results. 



XRef. solution (5x5 9ANS) 
▼5x5 4ANS, NLGEOM=1 
•5x5 4ANS, NLGEOM=2 



Transverse Deflection at Point b, mm 


(b) 4ANS results. 


Figure 15. Hinged Cylinder: Load vs. Free Edge Point Deflection 
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XRef. solution (5x5 9ANS) 
^5x5 AQ4 (ES24/A4S1) 
*7x7 AQ4 (ES24/A4S1) 



Transverse Deflection at Point a, mm 

(a) AQ4 results. 

X Ref. solution (5x5 9ANS) 
▼ 5x5 4ANS, NLGEOM=1 
• 5x5 4ANS, NLGEOM=2 


Transverse Deflection at Point a, mm 
(b) 4ANS results. 

Figure 16. Hinged Cylinder: Load vs. Center Point Deflection 
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Figure 17. Geometry of the Raasch challenge problem. 


Raasch Challenge Problem - The Raasch Challenge problem grew out of a presentation 
at the Structures Technical Forum at the 1 990 MSC World Users’ Conference. The geome- 
try is essentially that of a “hook” as shown in Figure 17. The hook is a thick curved strip 
clamped at one end (all degrees of freedom constrained to zero) and loaded by a unit in- 
plane shear load in the width (or z) direction at the other end. In this paper, the unit shear 
load is treated as a uniformly distributed shear load of 0.05 Ib/in. (total shear force of 1 .0 
lb). The hook has essentially two different curved segments that are tangent at their point 
of intersection. One segment spans an opening angle of 60-degrees and has a radius of 
14 inches. The second segment spans an opening angle of 150-degrees and has a radius 
of 46 inches. Both segments have a thickness t of 2.0 inches and a width w of 20 inches 
(aspect ratio w/t equal to 10 - limit for first-order, shear-deformation theory). The hook 
material is isotropic with an elastic modulus of 3300 psi and a Poisson’s ratio of 0.35. 

Early MSC/NASTRAN results by Harder (1 991 ) for this problem indicated that as the mesh 
of QUAD4 elements was refined, the solution did not converge unless the transverse shear 
flexibility (controlled by the parameter MID3) was suppressed. Recommendations to users 
at that time were to suppress the transverse shear flexibility (set MID3=0 on the PSHELL 
card) and also to suppress the drilling freedoms caused by the facet shell geometry approx- 
imations (used a large value for the artificial stiffness parameter K6ROT). The MSC ele- 
ment developers then proceeded to developed the QUADR element which included a dril- 
ling freedom at each node and also created the unique surface normal option for the 
low-order shell elements in order to resolve this behavior and to improve element perfor- 
mance of their low-order quadrilateral elements overthat reported by MacNeal and Harder 

(1985). 

At the 1995 MSC World Users’ Conference, Hoff et al. (1995) presented new results for this 
problem using both the QUAD4 and QUADR shell elements from MSC/NASTRAN Version 
68.2. The QUAD4 element is a 4-node flat, quadrilateral, displacement-based shell ele- 
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ment without any drilling freedoms. The QUADR element is related to the QUAD4 and in- 
cludes drilling freedoms at the nodes. Both shell elements can accommodate transverse 
shear flexibility. These new results also indicate that as the mesh is refined, the solution 
does not converge for either shell element unless the transverse shear flexibility is sup- 
pressed. The addition of drilling freedoms to the element formulation (i.e., the QUADR ele- 
ment) did provide some improvement in the element performance, but the solutions ob- 
tained using successively refined meshes still did not converge. The value of the tip 
deflection in the direction of the load obtained for the most refined mesh using the QUADR 
element was still more than twice their reported converged value (given as 5.01 2 inches 
by Hoff et al., 1995) and represents an improvement over the value obtained with the 
QUAD4 element which was nearly five time the converged value. Results were also pres- 
ented using a new surface normal option (SNORM parameter) which activates the creation 
of unique grid point normals for adjacent shell elements. Details of this option were not 
presented by Hoff et al. (1 995). With surface normals active and transverse shear flexibility 
suppressed, the solutions obtained without and with nodal drilling freedoms (QUAD4 and 
QUADR elements, respectively) are essentially the same and exhibit rapid convergence. 
In addition, the solutions obtained with surface normals active and also transverse shear 
flexibility included converge rapidly to their reference solution for both elements which cor- 
responds to a value approximately 5% more flexible than that obtained with transverse 
shear flexibility suppressed. Hoff et al. (1995) concluded that unique surface normals at 
grid points are needed to improve the element performance of the MSC low-order quadri- 
laterals. 

The results for various spatial discretizations and different shell elements shown in Figure 
18 are given in Table 3 wherein the isotropic properties of the material are used. That is, 
the transverse shear moduli are equal to the in-plane shear modulus (i.e., 
r - G - a = £ ) The reference value for the tip deflection used herein is 4.9352 

inches obtained using a refined 20x136x2 mesh of 8 node, assumed-stress hybrid solid 
brick elements (8_HYB). This value represents a 1 .6% stiffer solution than that given by 
Hoff et al. (1995). In all shell element models, the auto_dof_sup option is true, and the 
option auto_triad is false. The shell elements which neglect transverse shear flexibilities 
(4_STG, 3 DKT) converge to a solution that is slightly stifferthan the solution obtained us- 
ing a refined mesh of three-dimensional solid elements. This behavior indicates that trans- 
verse shear flexibilities can at most account for only 5% of the tip deflection. However, the 
solutions obtained with elements which do account for transverse shear flexibilities do not 
converge as the mesh is refined. The MIN3 and AQ4 elements give somewhat better solu- 
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(a) 1 x 9 mesh of elements 


(c) 5x34 mesh of elements 


(d) 10 x 68 mesh of elements 






(e) 20x136 mesh of elements <«> 20x136x2 mesh of elements 

Fiaure 18. Two-dimensional 4-node shell quadrilateral finite element models 
considered and three-dimensional 8-node brick finite element 
model used as reference solution. 
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tions than the 4_ANS and 4_HYB elements but still predict more than twice the tip deflec- 
tion of the three-dimensional solution for the most refined mesh. The MIN3 element sup- 
presses all drilling freedoms at the element level (artificial drilling stiffness), while the AQ4 
element includes the drilling freedoms to improve the in— plane displacement field approxi- 
mations. By including the drilling freedoms, the approximation for the in-plane displace- 
ment field component normal to an element edge is raised from linear to quadratic along 
that edge. Shell elements without transverse shear flexibility in COMET (4_STG and 
3 DKT) appear to converge to an appropriate value (i.e., 5% stiffer than the three-dimen- 
sional solution). Shell elements with transverse shear flexibility do not appear to converge. 
Surface normals or new computational nodal triads influence the results obtained using 
COMET but not as successfully as reported by Hoff et al. (1995). If the hook is made thin, 
then all shell elements tested converge 


Table 3. Curved Hook Results for In— plane Shear Loading 


Mesh of 
Elements 

Normalized Tip Deflection in Direction of Load 

Element Type 


Bh, XI» L 

4_ANS 

4_STG 

AQ4 

4_HYB 

MIN3 

3_DKT 

8_HEX 

8_HYB 

1x9 

0.9411 

0.9061 

1.1382 

1.1562 

1.0666 

0.8481 

0.5041 

0.6922 

3x17 

1.0335 

0.9398 

1.0375 

1.1678 

1.1341 

0.9323 

0.8136 

0.9432 

5x34 

1.2559 

0.9512 

1.1064 

1.3858 

1.2588 

0.9478 

0.9275 [ 

0.9738 

10x68 

2.0618 

0.9541 

1.3492 

2.1767 

1.6357 

0.9532 

0.9736 

0.9890 

20x136 

4.8111 

0.9548 

2.2529 

4.9045 

2.7760 

0.9479 

0.9932 

1.0000 


Results are normalized by solution obtained using a 20X136X2 mesh of 8 node, assumed— stress hybrid solid brick 
elements (8_HYB): 4.9352 for w/t=10. auto_dof_sup=true; auto_triad= false 
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Boeing HSCT Model - As a large-scale application study, a finite element model of the 
High Speed Civil Transport (HSCT) was analyzed using these elements. The analysis 
model and cases are company proprietary and specifics are not available. Ms. Tina Lotts 
of AS&M under contract for the Computational Mechanics Branch made these analysis 
runs and provided the normalized results shown in Table 4. The results are normalized by 
the results obtained by using MSC/NASTRAN as the finite element system. 


Table 4. HSCT Results for Various Load Cases 


Load Case 

COMET 

ES5/E410 

COMET 

ES24/A4S1 

1 

1.0692 

1.1183 

2 

1 .0333 

1.1142 

3 

0.9823 

1.0000 

4 

0.9755 

0.9895 

5 

1.0112 

1.0128 

7 

1 .0038 

1.0199 

8 

0.9830 

0.9830 

9 

1.0075 

1.0262 

10 

0.9823 

1.0000 
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RESEARCH PRODUCTS 

Under NASA Grants NAG-1-1 374 at Clemson University and NAG-1-1 505 at Old Domin- 
ion University, the following research products have been developed. The element re- 
search was facilitated by the use of the COMET software system and in particular the Ge- 
neric Element Processor or GEP. Two element processors were developed: ES22 for the 
2-node beam element and ES24 for the 4-node quadrilateral shell element. These pro- 
cessors have been transferred to NASA and are included in the baseline version of COM- 
ET. Three master’s and one doctoral students completed their degree requirements with 
full or partial support from this grant. Three conference papers and six journal papers have 
been published related to work sponsored by this grant with two additional papers still un- 
der review. Copies of each thesis or dissertation and each conference paper or journal ar- 
ticle have already been given to the technical monitor. 

Software 

• Processor ES22 - three-dimensional elastic beam element for linear stress, 
buckling, free vibration, and geometric nonlinear analysis using the low-order 
corotational formulation. 

• Processor ES24 - two-dimensional elastic quadrilateral shell element for linear 
stress, buckling, free vibration, and geometric nonlinear analysis using the low- 
order corotational formulation. 

Grant Review Presentations 

• May 28, 1992 

• June 30, 1994 

• January 6, 1995 

• December 12, 1995 

Graduate Students 

• Govind Rengarajan, “Assumed-Stress Hybrid Shell Elements with Drilling De- 
grees of Freedom,” Clemson University, August 1993, MS in Engineering Me- 
chanics. 

• V. R. Deshpande, “Development of an Assumed-Stress Hybrid C° Beam Ele- 
ment,” Clemson University, August 1993, MS in Engineering Mechanics. 

• David R. Oakley, “Concurrent Adaptive Dynamic Relaxation Algorithm for Non- 
linear Hyperelastic Structures,” Clemson University, May 1994, PhD in Engi- 
neering Mechanics. 

• W. Scott Carron, “Assumed-Stress Hybrid Beam Element for Nonlinear Elastic 
Stress Analysis,” Old Dominion University, May 1995, MS in Aerospace Engi- 
neering. 
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Conference Papers and Presentations 

• Oakley, D. R. and Knight, N. F., Jr. “Nonlinear Structural Response using Adap- 
tive Dynamic Relaxation on a Massively Parallel Processing System,” Proceed- 
ings of the Fifth Conference on Nonlinear Vibrations, Stability, and Dynamics of 
Structures, June 12-16, 1994, VPI&SU, Blacksburg, VA. 

• Carron, W. S., “Static Collapse of Circular Arches,” Presented at the 1 995 AIAA 
Mid-Atlantic Regional Student Conference, April 1995, NASA Langley Research 
Center (Second Place Winner based on written paper and oral presentation). 

• Knight, N. F., Jr., “The Raasch Challenge for Shell Elements,” AIAA Paper No. 
96-1369. Proceedings of the 37th AIAA/ASME/ASCE/AHS/ASC Structures, 
Structural Dynamics, and Materials Conference, Salt Lake City, UT, April 1 5-1 7 , 
1996, pp. 450-460. 

• Knight, N. F., Jr. and Carron, W. S., “Static Collapse of Elastic Circular Arches,” 
AIAA Paper No. 96-1 648. Proceedings of the 37th AIAA/ASME/ASCE/AHS/ 
ASC Structures, Structural Dynamics, and Materials Conference, Salt Lake City, 
UT, April 15-17,1 996, pp. 1 374-1 383. 


Journal Papers 

• Rengarajan, G., Aminpour, M. A., and Knight, N. F., Jr., “Improved Assumed- 
Stress Hybrid Shell Element with Drilling Degrees of Freedom for Linear Stress, 
Buckling, and Free Vibration Analyses,” International Journal for Numerical 
Methods in Engineering, Vol. 38, No. 11, 1995, pp. 1917-1943. 

• Rengarajan, G., Knight, Jr., N. F., and Aminpour, M. A., “Comparison of Symbolic 
and Numerical Integration Methods for an Assumed-Stress Hybrid Shell Ele- 
ment,” Communications in Numerical Methods in Engineering, Vol. 11 , No. 4, 
April 1995, pp. 307-316. 

• Oakley, D. R., and Knight, N. F., Jr., “Adaptive Dynamic Relaxation Algorithm for 
Nonlinear Hyperelastic Structures: Part I. Formulation,” Computer Methods in 
Applied Mechanics and Engineering, Vol. 126, 1995, pp. 67-89. 

• Oakley, D. R., and Knight, N. F., Jr., “Adaptive Dynamic Relaxation Algorithm for 
Nonlinear Hyperelastic Structures: Part II. Single-Processor Implementation,” 
Computer Methods in Applied Mechanics and Engineering, Vol. 126, 1995, pp. 
91-109. 

• Oakley, D. R., Knight, N. F., Jr., and Warner, D. D. “Adaptive Dynamic Relaxation 
Algorithm for Nonlinear Hyperelastic Structures: Part III. Parallel Implementa- 
tion,” Computer Methods in Applied Mechanics and Engineering, Vol. 126, 
1995, pp. 111-129. 

• Oakley, D. R., and Knight, N. F., Jr., “Nonlinear Structural Response using Adap- 
tive Dynamic Relaxation on a Massively-Parallel-Processing System,” Interna- 
tional Journal for Numerical Methods in Engineering, Vol. 39, No. 2, 1996, pp. 
235-259. 

• Knight, N. F., Jr., “The Raasch Challenge for Shell Elements,” Submitted to AIAA 
Journal, February 1996. 

• Knight, N. F., Jr. and Carron, W. S., “Static Collapse of Elastic Circular Arches,” 
Submitted to AIAA Journal, February 1 996. 
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CONCLUDING REMARKS 

Assumed-stress hybrid beam and shell elements have been formulated, implemented, 
and tested. The implementation was performed in COMET using the generic element pro- 
cessor. The beam element is available in element processor ES22 while the shell elements 
are available in element processor ES24. Numerical results for each element indicate that 
both element types are accurate and robust. Excellent correlation with published solutions 
and/or exact analytical solutions has been demonstrated for flat and curved structures, for 
thick and thin structures, and for linear stress, buckling and free vibration problems. Exten- 
sions to geometrically nonlinear problems has been achieved for both the beam and shell 
element by using the low-order corotational formulation. Excellent results have been ob- 
tained in comparison to other published solutions for nonlinear problems. However, as the 
deformation mode becomes more complex, more and more elements are required to mod- 
el the behavior accurately. For example, as the number of buckle halfwaves increases, the 
number of elements also increases due to the rapid change in the deformation mode. Aug- 
menting the linear strain-displacement relations with the nonlinear terms may relieve this 
requirement; however, the computational effort at the element level will increase signifi- 
cantly. 
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